﻿using System;

namespace iStudio.Library.Dsp
{
    public static class Convolution
    {
        /// <summary>
        /// Convolution using FFT
        /// </summary>
        /// <param name="signalOne"></param>
        /// <param name="signalTwo"></param>
        /// <returns></returns>
        public static float[] ConvolveFft (float[] signalOne, float[] signalTwo)
        {
            //
            // Make array's same length, to be able to perform FFT
            //

            if (signalOne.Length < signalTwo.Length)
            {
                Array.Resize(ref signalOne, signalTwo.Length);
            }
            else if (signalTwo.Length < signalOne.Length)
            {
                Array.Resize(ref signalTwo, signalOne.Length);
            }

            ComplexF[] signalOneComplex = new ComplexF[signalOne.Length];
            ComplexF[] signalTwoComplex = new ComplexF[signalTwo.Length];

            for (int i = 0; i < signalOne.Length; i++)
            {
                signalOneComplex[i].Re = signalOne[i];
                signalTwoComplex[i].Re = signalTwo[i];
            }

            Fourier.MyFft (true, signalOne.Length, signalOneComplex);
            Fourier.MyFft (true, signalTwo.Length, signalTwoComplex);

            //Fourier.Fft (signalOneComplex, signalOne.Length, FourierDirection.Forward);
            //Fourier.Fft (signalTwoComplex, signalOne.Length, FourierDirection.Forward);

            ComplexF[] signalFinalComplex = new ComplexF[signalOne.Length];

            for (int i = 0; i < signalOne.Length; i++)
            {
                signalFinalComplex[i] = signalOneComplex[i] * signalTwoComplex[i];
            }

            //Fourier.Fft (signalFinalComplex, signalFinalComplex.Length, FourierDirection.Backward);
        
            Fourier.MyFft (false, signalOne.Length, signalFinalComplex);

            float[] signalFinal = new float[signalFinalComplex.Length];

            for (int i = 0; i < signalFinal.Length; i++)
            {
                signalFinal[i] = signalFinalComplex[i].Re;
            }

            Effect.Normalize (signalFinal);

            return signalFinal;
        }

        /// <summary>
        /// A very simple mono convolution algorithm - NAudio Algorithm
        /// </summary>
        /// <remarks>
        /// This will be very slow
        /// </remarks>
        public static float[] Convolve(float[] input, float[] impulseResponse)
        {
            float[] output = new float[input.Length + impulseResponse.Length];

            for (int t = 0; t < output.Length; t++)
            {
                for (int n = 0; n < impulseResponse.Length; n++)
                {
                    if ((t >= n) && (t - n < input.Length))
                    {
                        output[t] += impulseResponse[n] * input[t - n];
                    }
                }
            }
            Effect.Normalize(output);
            return output;
        }
    }
}
